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ABSTRACT 


The objective of this thesis is to develop segmentation methods for multichannel 
and single channel images, and compare these methods. The segmentation algorithms 
are based on a linear model for the image textures and on inverse filtering to estimate 
the image textures and their regions. Two specific methods are compered 1) A 
multichannel filtering algorithm that simultaneously models the three separate signals 
representing the intensity of red, green, and blue as a function of spatial position and 
2) A single channel model applied to a combined image resulting from performing a 
Karhunen-Loeve transformation on the three signal components. Results of the 
multichannel image segmentation and the Karhunen-Loéve transformed one-channel 


image segmentation are presented and compared. 
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I. INTRODUCTION 


Segmentation techniques are among the most important considerations in the 
development of the automated image processing systems. Two related algorithms using 
2-D linear prediction models and the Karhunen-Loéve transformation for multichannel 
and color image segmentation are developed and compared in this thesis. 

The purpose of segmentation is to partition an image into a set of simpler 
homogeneous regions. The regions may consist of different gray level, different 
textures, colors, etc. Im some cases an ” image ” may consist of several spectral 
components. For example, a color image consists of three separate signals representing 
the intensity of red, green, and blue as a function of spatial position. If we represent 
each of these signals by functions F_ (n,m), F, (n,m), and C (n,m), the image is 


represented by a vector quantity 


e (n,m) 
F (n,m) = F, (n,m) (1.1) 
fe (n,m) 


where n and m represent spatial coordinates. We call such an image, consisting of 
more than one two dimensional (2-D) signal, a multichannel image. 

In this thesis, a method based upon linear prediction is evaluated experimentally. 
This method has been developed [Ref. 1] for monochrome images and extended to 
color images [Ref 2]. That method uses maximum likelihood (ML) and maximum a 
posteriori (MAP) estimation to segment multichannel images into regions of similar 
textures. The linear prediction is a filtering of the multichannel image to estimate the 
gray level at a particular spatial coordinate from the gray levels at neighboring 
positions. It is inplemented as a 2-D linear filtering operation. The algorithm uses a 
previously-determined set of parameters corresponding to the mean of the data in each 
channel, the covariance matrix of the prediction error, and the weighting coefficients of 
the estimation filter for each specific texture type. 

The method discussed above is compared to a variant of this method based on 


the Karhunen-Loeve (K-L) transformation. The K-L transformation allows the several 


components of an image to be combined into a single image that retains most of the 
energy in the original image. Hunt and Kubler [Ref. 3] found that for image 
restoration, Karhunen-Loeve transformation followed by single channel image 
processing worked nearly as well as multichannel image processing. It was desired to 
see if the Karhunen-Loéve transformation would be equally effective for segmentation. 
In this part of the work, the K-L transformation has been developed to reduce the 
3-channel color problem to a I-channel problem and the segmentation was performed 
for the one-channel image. Karhunen-Loéve transformation is based upon the 
Statistical characteristics of an image. The advantage of this approach is 
computational savings; only about one ninth the number of computations is required 
for this method. 

The remainder of this thesis is organized as follows. Chapter II discusses the 
model and the algorithms used to perform the multichannel image segmentation 
employing techniques of linear prediction [Ref: 4]. A general class of linear filtering 
models for texture is first presented. An algorithm is then developed to estimate the 
filter parameters from a multichannel image. Then, the multichannel image 
segmentation algorithm is described. 

Chapter III presents the models and the algorithms to perform the Karhunen- 
Loéve transformation and one-channel image segmentation. First, the algorithms for 
determining the eigenvectors and the eigenvalues of the correlation matrix are 
developed. Then, the transformation using a 3-channel image is presented. Finally, 
the one-channel image segmentation algorithm is discussed. The examples 
demonstrating the application of the segmentation methods for color images are 
presented and compared in Chapter IV. Chapter V has the conclusions about the 
multichannel image segmentation and one-channel image segmentation. 

In Appendix A, the Relaxation Method is described briefly as an alternative to 
the maximum a posteriori region estimation for monochrome images. Results are 
compared with the MAP method. In each of Appendices B and C, the description and 
use of the computer programs for the multichannel image segmentation and the one- 
channel image segmentation are presented. The computer program for multichannel 
image segmentation is contained in Appendix D. The Karhunen-Loéve transformation 
and one-channel image segmentation computer programs are contained in Appendix E. 
These programs are written in FORTRAN, compiled using Version 4.2 under the VAX 
/ VMS Version 4.1 operation system. 


Il. IMAGE SEGMENTATION USING A MULTICHANNEL FILTERING 
MODEL 

In this chapter, a multichannel image segmentation algorithm based upon a 2-D 
linear filtering model is presented. The multichannel images used in this work are color 
images with three channels representing the red, green, and blue components. The 
linear filtering model is used to develop approximate expressions for the multivariate 
probability density functions in terms of the filter error residuals for the entire set of 
points representing an image. The density expressions are used in the formulation and 
solution of the multichannel image segmentation problem. It is assumed that the 
multicharinel images contain multiple regions of homogeneous texture. This is found to 
be the case in dealing with aerial photographs of natural terrain. . 

' The problem of multichannel image segmentation is addressed as an estimation 
problem for two regions of texture. Maximum likelihood (ML) and maximum a 
posteriori (MAP) region estimates using a Markov random field to model region 
transitions are developed. 

This chapter consists of two sections. The first section describes the linear 
filtering model and develops an expression for the image probability density function in ~ 
terms of filter error residuals. The last section deals with the algorithm that is 
developed for the texture estimation and the multichannel image segmentation. The 


results of texture estimation and image segmentation are presented in Chapter IV. 


A. MODEL DEVELOPMENT 

In this section, a 2-D multichannel autoregressive-moving average (ARMA) model 
[Ref. 5] is first discussed. Then, we concentrate on the multichannel autoregressive 
(AR) model with Gaussian white noise inputs [Ref 6]. The development parallels that 
in [Ref. 1]. A multichannel image is represented by a vector signal F" (n,m) where 
(n,m) are spatial coordinates and the superscript h is an index representing the texture 
type. A 2-D multichannel image is shown in Figure 2.1 . A texture of type h is then 
modeled in general by a multichannel ARMA process defined by 


Panos 


ig? (nym) =) yy An. Fh (n-i,m-j) + Dy B". yy (n,m) (2.1) 
i=0j=0 B 
(i,j) # (0,0) 
Fo (n,m) = F" (n,m) + g® (2.2) 


for h = 0,1, n = l,....,N, m = I,....,M, where A. and B. are set of filter weighting 
coefficient matrices of size K by K. W" (n,m) are a set of independent identically 
distributed zero-mean random variables, G is a constant representing the mean value of 
the image, and B is finite-extent mask covering the filtered points. F" (n,m), W" (n,m), 
and on are vectors of size K, the number of channels in the image. The matrices An. are 
the kev parameters in the linear model. For a first quadrant filter with a P X Q region 
of support, there are PQ An. matrices with An = [, an identity matrix. 

For the auto-regressive (AR) or all-pole model we have B". os on and so that 
Equation 2.1 reduces to 


P-1t@al 

F" (nm) = — )) YA”, 2° (n-im-i) + ¥ (n,m) (2.3) 
i=0j=0 
(i,j) # (0,0) 


If the vectors f, w, and g represent an ordered set of the corresponding image points, 


then Equation 2.2 and Equation 2.3 are written in a matrix formulation as 
A(f- g)=W-Ajf, (2.4) 


where A and A, are matrices whose nonzero elements are derived from the terms Aj in 
Equation 2.1 and f) represents a set of boundary conditions with support outside of the 
regions. Since the terms W (n,m) are independent with probability density function 
(PDF) p,, One can solve Equation 2.4 for W and express the multivariate probability 


density function for the image conditioned on the boundary values as 


10 


m (n,m) 


F, (n,m) 


F, (n,m) 


F (n,m) 


F (n,m) = F, (nm) 


Peter 


Ech Der) 


2-D Multichannel Image Model. 


Channel K 


I 
Pri r, (f lf,)= Tat Pw (A(f-g) + Ag fig ) (2.5) 


= II p,, ( £ (n,m)) 
(nym) é€ R 


where the notation E(n,m) is used to represent the ordered components of the vector 


A( f-g) + Ap f,. If the boundary conditions, fy , are temporarily ignored, then 
E (n,m) = A (f — g) | (2.6) 


and the terms E& (n,m) in Equation 2.5 are computed from 


P-1 Q-1 
E" (am) = YY AS, EF (o-i.m-) | (2.7) 


i=Oj=0 
(1,J) # (0,0) 


The filter of Equation 2.7 which computes E h (n,m) from jou (n,m) is referred to as 
the ‘prediction error filter’. One can think of Equation 2.7 as producing an estimate or 
prediction F (n,m) of the image at point (n,m) and then forming an error E (n,m) as a 
difference F (n,m) — FP (n,m). This process is known as 2-D linear prediction and ts 


fundamental to performing multichannel image segmentation. 


B. ALGORITHM DEVELOPMENT 

In the multichannel segmentation problem, it is assumed that the image consists 
of multiple connected regions of known texture types, but that the region boundaries 
and the number of regions are not known. The segmentation of the image is treated as 
a supervised learning problem, since the regions are considered to consist of known 
texture types. In this section, the multichannel image segmentation algorithm for 
textured images is discussed. 


An overview of the method is as follows. Given a multichannel image of each 
texture, filter parameters are estimated by computing the covariance matrix from a set 
of data and solving the Normal equations corresponding to the model of Equation 2.3. 
In this case the ‘correlation method’ of linear prediction is used to compute the 
covariance matrix. The filter parameters are derived from a statistical analysis of the 
textured images, because the image model discussed in the previous section is based 
upon statistical properties. 

Once the filter parameters are known the filters are used to perform the 
segmentation. The filter weighting coefficients are used to calculate the prediction 
errors E® (n,m) of two textures (nm). Then, a maximum likelihood (ML) region 
estimate is developed using the prediction errors and the covariance matrices for image. 
The ML estimate is used as a basis to determine an approximate maximum a posteriori 
(MAP) region estimate. The MAP region estimation utilizes an underlying Markov 
structure for the region statistics to produce an accurate segmentation. 

1. Filter Parameter Estimation Method 

The prediction error filter is a finite-extent impulse response (FIR) or 
nonrecursive filter with selectable mask size and quarterplane region of support. It is 
always stable. The inverse of the filter used in Equation 2.3 1s an all-pole filter (the AR 
filter). The filter parameter estimation problem requires calculating Me pee and A? 
by statistical analysis of data in an estimation window containing the desired texture, 
where the quantity M" is a mean vector of the average gray level of the image in each 


of the channels, and ae 1S a cOVariance matrix of the multichannel white noise 


Eh) = E(w (n,m) (W* (n,m)? | (2.8) 


the term a is also refered to as the prediction error covariance matrix, since in a 
linear prediction problem it represents the covariance of the quantities E (n,m) defined 
meequation 2.7 . Since aoe is not in general a diagonal matrix, we see that the ‘white 
noise’ is uncorrelated within each channel but correlated between channels. 

The covariance of the multichannel white noise (the prediction error 
covariance) and the filter coefficients will be obtained by estimating the correlation 


function of the image and solving a set of Normal equations as discussed below. The 


is 


correlation function itself is estimated from data in a window containing a sample of 
the desired texture. Two estimation windows are depicted in Figure 2.2 for two 
different textures in the multichannel image. The reader should realize that the 
estimation windows do not have to come from the tmage to be segmented; they can be 


selected fron) any image containing the same type of texture. 






eae 


A__t 


Figure 2.2 Typical estimation windows for two textures. 






a. Mean Vector Estimation 
In order to model the multichannel image by Equations 2.2 and 2.3 the 
inean vector of the multichannel image has to be estimated. Knowing M and selecting 
the stationary estimation windows of two desired textures of F m? the Zero mean 2-D 
multichannel image, F, appearing in Equation 2.3 can be obtained by subtracting the 


mean vectors from the multichannel image. Thus, Equation 2.2 becomes 


E" (n,m) = E" | (n,m) - w EB) 


I4 


where M" corresponds to G" in Equation 2.2, and the term M" in Equation 2.9 are 


computed from 


M", 
m= |M>, | (2.10) 
M", 
where 
x, 
o>, > £", am) : (2.11) 


NM n=X,;m=yY, 


mM" is the mean estimate for the k"® channel of the h texture image. The limits Xi 
X5 ; ng ; a represents the edges of the window which is of size N’ by M’. Therefore 
N =X,-X,;+1, M =Y5- Y, + landO< X, <X, SN, 08 Y,<Y, 
S M ,and h represents the two textures of the multichannel image. All variables used 


in Equation 2.11 are depicted in Figure 2.3 . 
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Figure 2.3 Selection of the Size of the Windows for Mean Vector Estimation. 





b. Correlation Function Estimation 
The correlation function of the zero mean, 2-D multichannel signal has to 
be calculated in order to estimate the multichannel white noise covariance or prediction 
error covariance, Le and the filter weighting coefficients, A®. . The theoretical 2-D 


matrix correlation function for lag ” 1,j “ 1s given by 


RO (ij) = (RB (-i, jy)! (2.12) 
= E[F" (nm). ( F* (n-im-j))! ] 


and can be estimated from the multichannel signal by 


l ny m4 
R® (ij) = ————Y. YF" (n,m) ( £* (n-i, m-j)? (2.13) 
N M (ies 10\r; Ian feot: 


where R' (1,j) is a matrix of size K by K, and n, ,n,,m, , m, are defined by 

0 <n, = max(X, ,.X, + i) < n, = min(X,,X, + i) S$ N, and 

0S m, = max(Y,, Y, + j) < m, = min(Y,, Y, +j) S M. 
This matrix correlation function is used to form a larger block Toeplitz covariance 
matrix which is used to estimate the filter parameters. This 1s discussed next. 

c. Filter Coefficients and Prediction Error Covariance 

The prediction error filter weighting coefficients and the prediction error 
covariance must satisfy a set of linear equations known as the Normal Equations when 
the multichannel image is represented by the model in Equation 2.3. 


Normal equations corresponding to Equation 2.3 can be written as 
LR Seen (2.15) 


where R is the correlation matrix for the signal, A is an appropriately ordered matrix 
of the filter coefficients, and S contains a single non-zero block L, which is the 
prediction error covariance. The matrix R has three levels of partitioning and for anv 
rectangular region of support is block Teoplitz with block Toeplitz blocks. 

For a first quadrant filter with a P X Q region of support, The Normal 


equations that define Equation 2.15 have the specific form 
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R(O) R(-1) . 
R(1) R(0) 


| R(P-1) R(P-2). . 


where 


R(i,0) Ril) . 


R(i,1) R(i,0) 
R (i) 


LR(i,Q-1) R(i,Q-2) 
RE (.i) 


.R(-P+1) Al s? 

/R(-P+2){ | A! 0 

= (2.16) 
R(O) J LAP} 0 | 
. RG,-Q+ 1) 
- R(i,-Q +2) 
(OBIT) 
R(i,0) | 


and where R(i,j) is the matrix correlation function described in Equations .2.12 and 2.13 


The quantities A' and S° are defined by 


Ay 0 
. A ] 
A! = (2.18) 
AL Q-1 
and 
Ly 
0 
S = (2.20) 
0 


+ 


A 


with Apy= I and where A; and the partitions of S° are matrices of order K, the 
number of channels of image. 
2. Multichannel Segmentation Method 

An overview of the segmentation is as follows. By using a Gaussian 
probability density for the white noise in Equation 2.5 , one can develop explicit 
estimates for the density functions in terms of the prediction errors E (n,m) . From this 
one can form the conditions for ML and MAP estimation of the regions in the image. 
The theory leading to the estimates is explained below. A brief intuitive explanation of 
the process is given here. 

As mentioned earlier the prediction error filters in Equation 2.3 can be 
considered as predicting the intensity of a pixel in each channel from data in the region 
adjacent to the vector of pixels. The prediction errors E" (n,m) are the outputs of the 
filters. In the segmentation algorithm the prediction error is normalized in an 
expression involving the corresponding prediction error covariance. These normalized 
errors are compared in an appropriate formula to obtain the ML region estimate. 
When an area of texture is processed by a filter that is not matched to the texture, the 
normalized prediction error can be expected to be high. When the same area is 
processed by the filter that is matched to the texture, the prediction error can be 
expected to be low. 

The ‘maximum likelihood’ region estimate, ML(n,m), of texture class is 
achieved based on the prediction error covariance and the error estimates of the two 
desired textures for each pixel in the multichannel image. The ML(n,m) region 
estimate assigns pixels to texture types without regard to the assignments of the 
adjacent pixels. Then, the ‘maximum a posteriori’ region estimate, MAP(n,m), of 
texture class is achieved for a pixel and a desired number of adjacent pixels of two 
textures of the ML region estimation result. The MAP region estimation uses the 
Markov model that refers to the above description. The form of the Markov model 
and ML and MAP region estimates are presented in detail below. 

a. Maximum Likelihood Region Estimation 

It is supposed that a multichannel image has many regions, but that each 


region contains only one or another of two texture types. Given these regions, one can 
write the Equation 2.5 as 


PER) = Mpy (E (am) i 
(n,m) 


(2.21) 


ll 
a 
— 
e 
a 
pe, 
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where the p, is the probability density function for the white noise source of type h, 
within region R, ,and q is the number of regions. When the white noise term is 


Gaussian with density function (mean 0 and covariance L,, ) 


l l T FE 
y);=- - — Wt ly Dae 
Py (W) Nw? c 4 5 = W ) ( ) 
2) pe | 


then, taking minus twice the log of Equation 2.21 and applying Equation 2.22 , we 


obtain for an N by M pixel multichannel image 


SOMME IR, Ry j.---,R,) 


q 
=F dE, (amt Et, 7p! CEL, (my) + in|EpI +. 
(nm) eéR, 
+P (£4, (n.m)}? £4, 7 ! £4, (n,m)] + In|E4, |) 
(n,m) € Ra Seve In 2t 
| me 25) 
d ) ° e ‘ ° 
=> ¥ (ei, (am)? (ety | gi, (nm + nf, |) 
i=l (n,m)eéR, 


— NM In2n 


For maximum likelihood estimation, the number of regions q and the regions 
themselves are considered to be deterministic parameters of the density function. An 
ML estimate for these parameters is obtained by choosing values that maximize 
Equation 2.21 or, minimize Equation 2.23 . Since NMIn2z7 is constant value, the 
Equation 2.23 is minimized if every point (n,m) in the multichannel image to a region 
R, of type h, such that the term in brackets is minimum. Thus, one can write a ML 


region estimation for two textures as 


0 
ML, (n,m) 2 ML, (n,m) (2.24) 
1 


where 
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ML, (n,m) = [E* (n,m)? poh, 1 [B (n,m)]+ Info" (2.25) 


for h = 0, 1, where the number above or below the inequality indicates the region 
class to which the pixel (n,m) is assigned. When the class is 1, the ML(n,m) is assigned 
1 for a white pixel. Otherwise ML(n,m) is assigned 0 for a black pixel. Since the ML 
region estimate assigns pixels to black and white without regard to the assignments of 
adjacent pixels, this algorithm produces a number of false assignments and a somewhat 
‘Spotty result. 
b. Maximum A Posteriori Region Estimation 

The ‘maximum a posterior’ (MAP) region estimation utilizes the Markov 
model to describe the occurance of regions in the image. The combination of the linear 
filtering model with the Markov model results in an algorithm to achieve a MAP 
region estimation. For MAP estimation the regions are considered to be random 
quantities, and we maximize the probability for a given set of regions conditioned on 
our observation of the multichannel image. From Bayes rule, the a@ posteriori 


probability can be written as 


Pr[R, Ry seeeR, | 2] = a ae 


Since the denominator of Equation 2.26 does not depend-on the regions, maximum a 
posteriori (MAP) estimate for the regions can be obtained if the R, are chosen to 
maxinuze the numerator 


P(E |R) Ro vctaRa) Po (Rye Romer ee (2.27) 


af 


One can define the ” state ” s(n,m) of point (n,m) as the region type to 
which that point has been assigned. In our development, the number of region types is 
assumed to be 2. Since the set of all possible state assignments for points in the image 
iS One-to-one with the set of all possible divisions of the multichannel image into 
regions, the region estimation problem can be viewed as one of estimating the states of 
the points. It is assumed that the state of a point is stochastically dependent on some 
adjacent set of states Sn, m i a symmetric support region, as shown in Figure 2.4 . 


Since S represents a chosen set of state assignments for all points in the multichannel 


image, Pr [ S ] denotes the joint probability that the points in the image take on a 
chosen set of state assignments. The support set defines a neighborhood structure L.e. 
all elements in the support set are neighbors of each other. It can be shown that if the 
set of states S is a Markov random field, then the probability of S can be factored as a 
product of terms of functions depending only on the “cliques” of the support setS, _ . 
The cliques are defined as groups of points such that each set of points are neienbor 
of each other according to the support set. For a Markov random field, the probability 
of S can be written as 


? 





Figure 2.4 State support region of the point s for MAP estimation. 


Po] S)] = II Pr[ s(n,m) | S ] (2:28) 
(n,m) Nn, m 


where the terms in the product are additive functions defined on the cliques of S| 
and the product is over all cliques in S. One simple acceptable form for the terms in 


the product is 


Pr [ s(n,m) | Sn, | ~ exp [s(nm) {a + B, (ttt) + B, (vt+v’) 
+ ¥, (utu’) + ¥, (wt+w’)}] (2.29) 


where t = s(n-l,m), t = s( n+1, m), Sm is the set of states as shown in Figure 2.4 , 
and D is a normalizing constant. One particular selection of the parameters, namelv @ 
a P= B, = ¥,; = Y> = 1, leads to a particularly simple algorithm. In this 
case, we have 
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| 
Pr 1S, 4, | = ——= expiey, (sGiee eal ace) (2.30) 
D (1) ES. om 


and 


1 


ae (2.31) 


de a 


The second term in the numerator of Equation 2.26 can be replaced by 


Pr[S] . Thus maximizing Equation 2.26 is equivalent to minimizing 
— In p( FR, peorRy ) = 2in Pai (2.32) 


One can define the MAP region estimate by combining Equations 2.22 , 
2: 269n 2-287, and) 2.3caas 


e 


rz (n,m) ]! (E17! [E! (n,m) 
(Hym) 


0 
] aa 
+ In| 2, | 2InPr{[l]S, als 
l 
ee (n,m)}? (£9 7 | [£° (n,m) | (2.33) 
0 = 
+ In) 22 ee nee Oe 


Por er aene Si m | in the form of Equation 2.30 and Equation 2.31 computing the 
COrimSaamee itr [hls m | is equivalent to counting the number of pixels in o m that 
have value ‘h’ and dividing by the total number of pixels in S a and multiplving by 
an appropriate normalizing factor, Ks! A larger state support region Sa, m iS depicted 
in Figure 2.5 as an example. The side of the S n, m Must be an odd number . Although 
it does not necessarily lead to the true MAP estimate we find it convenient in practice 
to maximize Equation 2.33 term by term. That is, we require to use Equation 2.33 


without sum. Equation 2.33 can be solved by iteration using the maximum likelihood 





lThe normalizing factor results because the pa a,f 


Equation 2.29 can bé scaled arbitrarily and still result in a Naeithdte aaa 
function. 
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state estimates obtained from Equation 2.24 as an initial set of states for MAP region 


estimate. 





Figure 2.5 A Set of States I’s, 0’s adjacent to S for MAP Region Estimate. 


The computational requirements of the MAP region estimation are reduced by storing 
the differences, 


MLD (n,m) = ML, (nm) — MLo (n,m) (2.34) 


incurred during the calculation of the ML region estimate. Substituting Equation 2.34 
in Equation 2.33 gives 


m 


a. 1 
MLD(n,m) — 2 In Pr{ |S, J+ 2inPr{Ols, ,1< 0 (2.35) 


0 


At each iteration terms Pr[h|S___ J are evaluated based on the values of the states 
9 
at the previous iteration. For our particular method of selecting Pr [{ h | S 


_ Equation 2.35 can be expressed as 


na 


M LD(n,m) 


2 (number of state 1 pixels) — ( number of pixelsinS, —) + | 
= KS 


SN Nae 
> 


number of pixels in S. i 
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(2.36) 


There are two important points in Equation 2.36 in order to perform the maximum a 
posteriori image segmentation accurately. The value of the convergence factor, KS, 
must be assigned properly. If KS is assigned too small, the segmentation may not 
remove improperly classified pixels. On the other hand, if KS is assigned too large, 
correctly classified pixels could be changed. In addition, the size of S_ m must be large 
enough. Otherwise, false assignments produced by the initial ML segmentation may 


not be removed. 


Ill. KARHUNEN-LOEVE TRANSFORMATION AND ONE-CHANNEL 
IMAGE SEGMENTATIO 

In this chapter, the models and relevant algorithms are presented to perform the 
Karhunen-Loéve (K-L) transformation [Ref. 3] and one-channel image segmentation 
utilizing the techniques of linear prediction [Ref. 1]. Since linear prediction techniques 
are presented in detail for multichannel image segmentation in the previous chapter, 
~ this discussion concentrates on the K-L transformation. The K-L model and algorithm 
are first developed to reduce the three-channel color problem to a one-channel 
problem. Then, a one-channel segmentation procedure is presented that is based on the 
Same model previously discussed. The results of the K-L transformed one-channel 
image segmentation are presented and compared with the multichannel image 


segmentation in Chapter IV. 


A. MODEL DEVELOPMENT 

The K-L transformation developed in this section is based on the statistical 
properties of an image. This transformation provides an energy compaction between 
channels of a color image. That is, most of the color image energy is compacted into 
one channel, and the transformed image channels are uncorrelated. If the multichannel 
image and transformed multichannel image are expressed in vector form, The K-L 


transformation is given by [Ref. 7] 
no) A (n,n) (3.1) 
where F ( n,m ) is the original multichannel image, Q ( n , m ) is the transformed 
multichannel image, and A is the K-L transformation matrix, whose rows are 
eigenvectors of the between-channel correlation matrix R defined by 
R= E[F(n,m)£!(n,m)}, (3.2) 


The between-channel correlation matrix of the transformed image is 


A= E[@(n,m)@! (n,m)} (3.3) 
=(AJ[R][AT! 


2 


where the matrix A is a diagonal matrix of the form 


yO) ae 
A= Ch nO (3.4) 
OCS 


and the A, represent eigenvalues of the between-channel correlation matrix. The 


eigenvalues are ordered such that 
A, 21,21, 20 (3.5) 


The importance of this property is that each eigenvalue hy is equal to the variance of 
the k ™ channel of the transformed multichannel image whose channels - are 
uncorrelated. Then, it is a well known property of multivariate statistics that the total 


variability of the color image has the form [Ref. 3] 
m= DA, (3.6) 


which relates the total variability to the decorrelated component variations, A, - One 
can observe that often the Ae values have a wide range of magnitudes, and the first 
component, A, , will be sufficient to approximate A7 with only a small percentage of 
error. This becomes the key idea for the use of the K-L transformation. Table 1 shows 
the energy distribution between the transformed color image channels of two test 
images. 

Indeed, in a 3-channel image one can often find that the first channel of the K-L 
transformed image is sufficient to account for 99 percent of all the variability. That is, 
it is typical to find that 


h, = 0.99 A7 (3.7) 


The Karhunen-Loéve transformation provides the best energy compaction [Ref. 8] and 


the advantage of this transformation is computational savings. We will see later that 


TABLE I 
ENERGY DISTRIBUTION BETWEEN TRANSFORMED IMAGE 
CHANNELS 


Percentage of Energy in Channels 
3 
99.13 0.84 0.03 
99.15 0.78 0.07 





one-channel image segmentation achieved bv processing only the first channel of the 
K-L transformed color image will be very close to the result obtained on the original 
image with a multichannel algorithm, but will be obtained at only one ninth of the 


computational cost. 


B. ALGORITHM DEVELOPMENT 
In this section, the algorithm to perform K-L transformation of color images is 
presented. Since the K-L transformation is based upon the between-channel correlation 
matrix, R , of the color image , the between-channel correlation matrix is first 
determined. Then the eigenvectors of the between-channel correlation matrix are 
computed. Finally, the K-L transformation is formulated using the transpose of the 
eigenvector matrix. 
I. Correlation Function Estimation 
In order to determine the K-L transformation of Equation 3.1 , the between- 
channel correlation matrix of the color image must first be estimated. Our estimate for 


the between-channel correlation matrix is given by 





1 N-I N-1 
R = YYE(n,m)(E(n,m))t (3.8) 
N? n=0 m=0 


oe 


where the image size is N by N, and the between-channel correlation matrix size is K 
by K. 
2. Karhunen-Loeve Transformation Matrix 
Since the rows of the K-L transformation matrix are the eigenvectors, E, of 
the between-channel correlation matrix, the eigenvectors must be calculated from the 
between-channel correlation matrix of the color image. The K-L transformation 
matrix is then obtained using the transpose of the eigenvector matrix. That is 


<—el —_ 


I 
[A ]= <7, a” (3.9) 


A 


<—¢!,—» 


3. Karhunen-Loeve Transformation 
Since the K-L transformation matrix satisfies Equation 3.3 , then the K-L 


transformation is represented by Equation 3.1 with A is given by Equation 3.9, where 


T 


e",,1 = 1,2, 3 are the eigenvectors of R. More explicitly, the components of @ are 


determined from the components of F at any pixel( n,m) as 


Q, (n,m) «—_!  —> Fy (nom) 
Q,(n,m) |=|<—2! ,»|IF, (n,m) (3.10) 
Q,(n,m) <=_—  ¢ 7a F3 (on gaas) 


C. ONE-CHANNEL IMAGE SEGMENTATION 

The one-channel segmentation procedures presented in this section are based on 
the same model that was described in detail for the multichannel image in Chapter II. 
A single channel image is used instead of a three channel image. That is, K is always 
equal to | in the Equations of Chapter II. As a result the vector and matrix quantities 
become scalars and the equations are simplified. 

In order to model the single image by Equations 2.2 and 2.3, the mean of the 
image is estimated using the Equation 2.11 . Then the correlation function is estimated 
in the same manner as in Section II-B.b where R(i.j) is a scalar instead of a matrix. 
This procedure is followed by estimation of the filter coefficients and the prediction 


error covariance. The equations in Section II-B.c are then used to determine the filter 
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coefficients, A. , and the the prediction error covariance, & w » now a scalar value. 
Finally, the one-channel segmentation algorithm 1s applied. The method is the same as 
that described in detail for the color images in Section II-B.2. However, instead of 
Equation 2.24 and 2.33 the following simplified relations are used for the maximum 
likelihood and the maximum a posteriori region estimates. A maximum likelihood 


region estimate for a single channel image of the pixel is given by [Ref. 1] 


S (£9 (nm)? 
t ae sane) “<4 ———sT 


WwW WwW 


( E! (n,m)? 


a cael nares) (3.11) 


—" 


and the maximum a posteriori region estimate is given by 


2 ; 
( E! (n,m) 
=a ei Danvers (AlN Sie] 





0 Ww 
Sem)? _, (3.12) 
< le eee nee OS | 


where E! and E® are the result of applying the linear predictive filters to the first 
channel of the K-L test image. 


IV. RESULTS AND COMPARISON OF THE METHODS 


In this chapter, the results of the multichannel image segmentation, the 
Karhunen-Loéve transformed one-channel image segmentation, and the K-L 
transformed 3-channel image segmentation are presented and compared. 

The digitized image size used in this work is 128 by 128 pixels with gray levels 
represented on a scale of 0 to 255 (8 bits). A digitized color photograph of a rural area 
containing trees (the green region) and fields (the yellow region) are shown in Figure 
4.1. A quarter-plane filter for each texture class (2 X 2 pixels) was designed and 
applied to the color image to achieve the multichannel image segmentation. A state 
support region of 7 by 7 pixels was used for MAP region estimation. The results of the 
maximum likelihood (ML) and the maximum a posteriori (MAP) segmentations of 
Figure 4.1 image are shown in Figures 4.2 and 4.3 respectively. The segmentation 
results show the field regions as black and tree regions coded as white. The ML result 
(Fig 4.2) is spotty, but the true tree and field regions are distinguishable. The MAP 
segmentation result (Fig 4.3) of Figure 4.1 image is quite clear. The MAP 
segmentation was not able to remove a few improperly classified points in the left side 
of the field region. Since there is an inherent ambiguity in the estimation of the region 
boundaries due to the finite size of the masks, the edges are expectedly somewhat 
rough. Figure 4.4 shows another color image of a rural area containing treés and 
fields. Figures 4.5 and 4.6 show the results of the ML and the MAP segmentation of 
the Figure 4.4 image. The ML estimation was achieved using the filters designed for 
the image of Figure 4.1 . The result of the ML segmentation (Fig 4.5) is again quite 
spotty, ie. Although two regions are discernible. The MAP algorithms segmented the 
regions quite accurately as shown in Figure 4.6 . The MAP estimates presented above 
converged after 10 iterations and KS was assigned to 100. 

In the results presented above the ML procedure produced a poor result with a 
lot of false regions. This is due to the lack of prior information about region 
connectivity. On the other hand, MAP estimation using the Markov model to represent 
region transitions produced results that was quite accurate. 

Figure 4.7 shows the first channel of the Karhunen-Loéve (K-L) transformed 


image of the Figure 4.1 . The image size is 128 by 128 pixels with the scaled gray levels 
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Within the intensity range of the display. The result of the one-channel ML 
segmentation of Figure 4.7 is depicted in Figure 4.8 . The one-channel ML 
segmentation result is quite spotty, but the two regions are perceptually discernible. On 
the other hand, the MAP segmentation result (Fig 4.9) of Figure 4.7 is quite smooth, 
although there are a few incorrectly classified points in both regions. 

The first channel of the Karhunen-Loéve transformed image of Figure 4.4 is 
shown in Figure 4.10 , and Figures 4.11 and 4.12 show the results of the one-channel 
ML and MAP segmentations of Figure 4.10 respectively. The one-channel ML 
segmentation result (Fig 4.11) has a lot of spots, but both segmented regions are 
distinguishable. The maximum a posteriori segmentation (Fig 4.12) of Figure 4.10 is 
quite smooth, but again there are a few incorrectly classified sets of pixels in both 
regions. , 

The results of multichannel image segmentation and the Karhunen-Loéve 
transformed one-channel image segmentations were presented consecutively. These 
results show that the ML estimation of the Karhunen-Loéve transformed single 
channel image is much more spotty than the ML estimation of the multichannel image. 
But, the MAP estimation results of the K-L transformed single channel images are as 
clear as the MAP segmentation results of the multichannel images. 

In summary the results presented in this chapter show that the best segmentation 
result is provided by the multichannel image segmentation method. However the 
segmentation results of the Karhunen-Loeve transformed single channel images are 
very close to multichannel image segmentation results, i.e. Because most of the color 


image energy (99 percent) is compacted into the first channel. 
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Figure 4.1 Two-Texture Color Image. 








Figure 4.2 ML Region Estimation of Figure 4.1 Image. 





Figure 4.3 MAP Region Estimation of Figure 4.1 Image. 
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Figure 4.4 Color Image Containing Two- Texture. 
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Figure 4.5 ML Region Estimation of Figure 4.4 [mage. 





Figure 4.6 MAP Region Estimation of Figure 4.4 Image. 
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Figure 4.7 First Channel of K-L Transformation of Figure 4.1 Image. 





Figure 4.9 MAP Region Estimation of Figure 4.7. 
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Figure 4.10 K-L Transformed Single Channel [mage of Figure 4.4. 
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Figure 4.11 ML Segmentation of Figure 4.10. 





Figure 4.12 MAP Segmentation of Figure 4.10. 
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V. CONCLUSIONS 


The segmentation of terrain images is an important part of image analysis 
methods for military and civilian applications. The work in this thesis utilized a 2-D 
stochastic linear filtering model and compared algorithms for multichannel image 
segmentation of color images. Two levels of structure were used for multichannel 
segmentation development. The fundamental structure based on the linear filtering 
concepts represents the texture in local regions of terrain. Superimposed on this 
structure is a Markov random field that describes transitions from one region type to 
another. The segmentation was considered as a region estimation problem and 
maximum likelihood and maximum a poSteriori region estimatation methods were 
developed. The ML region estimation produced a spotty result, but the MAP region 
estimation produced quite accurate results for the multichannel and single channel 
image. 

The other piece of work developed in this thesis was the Karhunen-Loéve 
transformation model that based on the statistical characteristics of color image. The 
one-channel image segmentation was then applied to the first channel of the 
Karhunen-Loeve transformed color image to see the effectiveness of the K-L 
transformation for segmentation. 

We observed that multichannel image segmentation results were quite accurate. 
Similarly the results of the K-L transformed one-channel image segmentation were very 
smooth. In summary the results of both segmentation methods were very close to each 


other, and the K-L transformation is very effective for segmentation. 
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APPENDIX A 
RELAXATION METHOD 


The Relaxation Method algorithm utilizes a set of iterative numerical techniques 
to compute a posterior probability of the pixel (n,m) from a prior probability of the 
same pixel and a set of prior probabilities of the adjacent pixels. The prior probabilities 
are estimated from a 2-D image data using the linear filtering model (see Equation 
ue 5). 

The Relaxation formula is defined [Ref. 9] by 


2.5 P..* (n, 
P.. K +1 (nm) = Avg ( es (A.1) 


$ k 
a P.. (n,m) 


where k is the number of the iteration, hie is the relaxation factor, P.. (n,m) are a set 
of prior probabilities, T is the number of textures, and S is the number of pixels. The 
updated estimates Ee +I (n,m) are obtained by averaging all of the terms in 


parentheses. The relaxation factor, hie , in Equation A.1 is given by 


ae 
hi? = E clii | sit) Pee (nm) (A.2) 


where c(i,j|s,t) is a nonnegative compatibility function whose value is small if the 
neighboring pixel is black when the estimated pixel is white, otherwise its value will be 
large. hie is used to update the probability Ps (n,m). Note that hie is large if the 
compatibilities c(i,j]s,t) are large and the probabilities Pe are high, otherwise hie will 
be small. 

The results of the Relaxation Method segmentation are shown in Figures A.1 
and A.2 . The Figure A.l presents the segmentation of Figure 4.7 with the 
compatibility c(i,x|s,x) is equal to 0.1, where x can be 0 or 1. The result is very spotty 
but both regions are discernable. The Figure A.2 gives the segmentation with the 
compatibility c(i,x|s,x) is equal to 0.9. This result is better than the previous result, but 


there are still several spots especially in the field region. Both results are obtained after 
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[0 iterations. Figure 4.9 shows the result of MAP segmentation of Figure 4.7. The 


MAP segmentation result is much more accurate than the Relaxation method 


segmentation. 





Figure A.1 Relaxation Method Segmentation with c = 0.1. 
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- 


Figure A.2 Relaxation Method Result with c = 0.9. 
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APPENDIX B 
FILTER PARAMETER_ESTIMATION AND MULTICHANNEL 
SEGMENTATION 
The interactive program, FLTR1I, estimates two sets of filter parameters from a 
2-D three-channel image. The corresponding algorithms were presented in Section 
II.B.1. In order to run this program, the user must input the required program 
parameters in the order listed below 
* The number of filter rows, P. Maximum is 4. 
* The number of filter columns, Q. Maximum is 4: 
* The number of rows, N, in each image channel. Maximum is 128. 
e number of columns, M, in each image channel. Maximum is 128. 
* The number of channels, K, in image. 
* The filenames of the image channels. 
* ‘The coordinates of the estimation windows of textures. 
* Two output filenames for the estimated filter parameters. 
The mean vectors of the data in the two estimation windows is first computed by 
FLTRI. Then, the program subtracts the mean from the image and determines the 
correlation functions. After calculating the correlation matrix, the program determines 


the covariance matrix, L , using the equation 


é 


[RIEBI = [Sy] | (B.1) 


° 


where B is a dummy matrix that has the same dimension with A matrix, Bog = { xy yi 
and Sy is all zero except for an identity matrix in its first partition. The covariance 


matrix must satisfy the relation 


Bale (B.2) 


WwW 


Finally, the filter weighting coefficients are estimated using 


[A]=[B].[=] (B.3) 
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MULTICHANNEL IMAGE SEGMENTATION PROGRAM 


The interactive program, SGMTI, segments a color image with two textures 
when given the three-channel image, the image dimension, and two sets of filter 
parameters. Again, the user must input the required program parameters to run the 

rogram. These parameters have to be in the order listed below : 
prog P 

* The number of rows, N, in the image channels. 

* The number of columns, M, in the image channels. 

* The number of filter rows, P. 

* The number of filter columns, Q. | 

felhe number of textures, I, in the image. 

* The number of channels, K, in the image. 

* The filenames of the image channels. 

* The filenames of the filter parameters. . 

* The output filename of the ML segmentation. 

* The output filename of the MAP segmentation. 

This program estimates the errors of two textures using Equation 2.7, then performs 
the ML segmentation and the MAP segmentation using Equations-2.25 and 2.33 . The 
convergence factor KS, and the size of S| must be assigned properly by user to 


perform the MAP segmentation accurately. 
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APPENDIX C 
KARHUNEN-LOEVE TRANSFORMATION 


The interactive program, KRLV, implements the Karhunen-Loéve transformation 
algorithms described in Section III.B. This program requires the user to assign the 
program parameters in the order listed below : 

* ‘The number of channels, K, in the image. 

* The number of rows and columns, N, in each channel. 

* The filenames of the image channels. 

* The output filenames of the transformed color channels. 

The correlation matrix is first calculated. Then the program calls the IMSL 
subroutine EIGRS to calculate the eigenvalues and the eigenvectors of the correlation 
matrix. Finally, the transformed color image is obtained by multiplying the transpose 
of the eigenvectors matrix by the original color image. The program scales the 


transformed color image for display on the COMTAL image prossessing system. 


ONE-CHANNEL IMAGE SEGMENTATION PROGRAM 


The program, FLTR2, calculates two sets of filter parameters from a 2-D single 
channel image. The user must input the required program parameters to run the 
program. These parameters except the number of channel, K, are given in Appendix B. 

The program, SGMT2, implements a single channel image segmentation. The 
required program parameters. given in Appendix B have to be assigned to run the 
program. Equations 3.11 and 3.12 are used to perform the ML region estimation and 


the MAP region estimation. 
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QAM aA aA 


APPENDIX D 


OMPUTER PROGRAMS FOR MULTICHANNEL IMAGE 
ee ae SEGMENTATION 


KRRRKKKRKKRRKRRRKRRERRRRRRRREKRRRKRRERRRRRRRRRRERERRRERRRRKRRERKRERRRRERARRRRERE 


; PROGRAM FLTR1 
GeeerURPOSE To develop two sets of filter parameters from a 2-D 
a - three-channel input image. These parameters are the 
* mean vector, the covariance, and the filter coeffici- 
: ents of the image. 

: REQUIRED IMSL ROUTINES 

- LEQT2F, LUDATN, LUELMN, LUREFN 

*k 

x 

x 

x 


IMPLEMENTED BY LTJG TIMUR KUPELI Nov 1986 


AO OO 


RERKEKRKKKRARREKRRKRKRKRRRKRRRKRKRRKARRRRRKRKRRKRRRRRKRRERKRKRKRRRRRRRKRRRRRRRRRRRRRRARKRERR 


AkkX VARIABLE DEFINITIONS ARRK 


BINPUT = Input |] Image data in byte format. 

FINPUT = { Input | image data in real*4 format. 

P , g = Rows, Columns of the linear filter.: 

N , = Rows, Columns of the image data. 

K = The number ‘of channels in image. | 

z -TUlemiuuDerwon Lilters for processing. 

IMAGE = : poeue }_Filename of the image channels. 
FILTER = [ Output ] Filename of the filter parameters. 
MEAN = Array of mean vectors of estimation windows. 
R = The correlation matrix. 

IMATRX = Identity matrix. - 

KW = The covariance matrix 

A =Siemeumnter coefricients matrix. 

ISIZE = The estimation window size. 

NO, MO = Row, column delay(shift) of correlation. 


This program uses 2 by 2 pixels filter. If user wants to use 
gLEterent Sige f£iiter, the dimensions of the following 
variables must be modified. 

SI, A, WAREA 

e, for 4 Tas cae filter, the dimensions must be 


R 
BOr exemp 1 ; 
(48,48), SI ), A(48,3), WKAREA(2448) 


INTEGER*4 P,Q,ROW,COL,STARTN(2),STARTM(2) ,ENDN(2),ENDM(2), 


* RNROW, RNCOL, RROW,CCOL,X,Y,K,1IER,QK,PKO,L,J,JJ,JJd, 
x HNNO ,HMMO,HCOL,HROW, LNNO, LNMO, LCOL, LROW,RN,RM,T, 
* NO,MO,RRN,RRM,ROW1,COL1,MMO,NNO,IDGT,I 

REAL*4 ey ee aS AREA (180) . 
* R(12,12),$1(12,3) , IMATRX(3,3),B00(3,3),SUM1,SUM2,SUM3, 
* KW(3,3),A(12,3),1SIZE 


CHARACTER*50 IMAGE(3),FILTER(2) 
BYTE BINEUT L288) 
T = 2 

GET PROGRAM INPUT PARAMETERS. 
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LYFE at 

10 FORMAT(' ENTER NUMBER OF FILTER ROWS FROM 2-4 DESIRED :',§$) 
READ 11,P 

it FORMAT (13) 


TYPE 


12 FORMAT” ENTER NUMBER OF FILTER COLUMNS FROM 2-4 DESIRED :',$) 
READ 11,0 


i ees 

i3 FORMAT(' ENTER NUMBER OF ROWS IN IMAGE ',$) 
READ 11,N 
TYEE 


14 FORMAT (' ENTER NUMBER OF COLUMNS IN IMAGE ',$) 
READ 11,M 


15 FORMAT(* ENTER NUMBER OF CHANNELS [ MAX = 3 ] IN IMAGE ;',$) 
READ(%*,11) K 
GET THE MULTICHANNEL IMAGE 
po1l6éJ=1,K 
WRITE(*,17) J 
19 FORMAT(' ENTER FILENAME OF IMAGE CHANNEL ' I3,$) 
READ (*, 18) IMAGE (J) 
18 FORMAT (A5O) 
CONVERT THE IMAGE FROM BYTE FORMAT. TO REAL*4 FORMAT 
OPEN(UNIT=1 , FILE=IMAGE(J) ,STATUS='OLD',ACCESS = 'DIRECT') 


DO 180 ROW = 
READ (2+ ROW) SET ae ,COL=1 ,M) 


DO 181 
TEMP = ‘BINPUT (COL) 
IF (TEMP. LT .0.0) @rReEn 
TEMP = TEMES 2oG 


END IF 
FINPUT(ROW,COL,J) = TEMP 


Qa AaAAaANQ 


OY AAN 


181 CONTINUE 
180 CONTINUE 
CLOSE (UNIT=1) 


16 CONTINUE 


GET 'T' AREAS FOR WHICH FILTERS ARE DESIRED AND OUTPUT FILENAMES 
FOR EACH AREA'S FILTER COEFFICIENTS AND COVARIANCE MATRIX. 


DO 19 ~«+XI a 

WRITE (*, 30) I 

20 FORMAT(' ENTER UPPER-LEFT ROW FOR AREA ',12,':',$) 
READ(*,11) STARTN(I) 


WRITE(*,21) I 
Ze FORMAT ( ' ENTER UPPER-LEFT COLUMN FOR AREA! Ae ae) 
READ(*,11) STARTM(I) 


WRETE (4 522) 1 
22 FORMAT ( ' ENTER LOWER-RIGHT ROW FOR AREA! eRe es) 
READ(*,11) ENDN(I) 


WRITE(*,23) I 
23 FORMAT(' ENTER LOWER-RIGTH COLUMN FOR AREA',I2,':',$) 
READ(*,11) ENDM(I) 


AAA 4a 
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WRITE(*,24) I 
24 FORMAT (! ENTER OUTPUT FILE-SPEC FOR MWiiteR ~ 12,53, 9) 
READ (*,18) FILTER(I) 


FIND THE MEAN VECTOR OF ESTIMATION WINDOW , T, AREAS OF IMAGE 
CHANNELS. THE MEAN VECTOR CONSISTS OF THE MEAN FOUND EACH CHANNEL 


aaa See - STARTN(I) + 1)*(ENDM(I) - STARTM(I) + 1) 
UM ° 


ANNAN 


S 

Nr 
1 
Oo 

Oo 


po aa ENDN(I) 
ARTM(I), ENDM(I) 
SUM + FINPUT eel 
SUM2Z + FINPUT et 
SUMGes FINPUT(L,J,3 


ng 


0 
DO 25 L = ST 
DO 2 


26 CONTINUE 
25 CONTINUE 
MEAN I, 1) = SUM1 / ISIZE 
MEAN 2 = SUM2 / ISIZE 
MEAN a = SUM3 / ISIZE 
WRITE (*, ,27)_ (MBAN(I, IT) 
27 ¥FORMAT('’ MEAN:',3F9.3,$) 


CORRECT THE IMAGE TO BE ZERO MEAN 


Dowel 
DO 28 L= STARTN(T) , ENDN(I) 
DO 29. LL = STARTM(I),ENDM(I) 
FINPUT(L,LL,J) = FINPUT(L,LL,J) - MEAN(I,J) 
29 CONTINUE 
28 CONTINUE 
271 CONTINUE 


DETERMINE THE 2-D CORRELATION FUNCTION OF THE IMAGE CHANNEL. 
Ee oN MATRIX APPROPRIATE TO THE INPUT INPUT PARAMETERS 
, K, Q 


WRITE(*,291) I . 

291 FORMAT (! CORRELATION.MATRIX'! ,12,$) 
OES sence ~ 0 
OK=QO%*K 

DO 30 " RNROW = 1,P 
= OK * (RNROW-1) 

50 31 RNCOL = 1,P 
No_= RNROW - RNCOL 


DO "33 RMROW = 1,0 
RN = K & (RMROW -1)+X 
DO 33 RMCOL = 1,0 
MO RMROW - “RMCOL 
RM * (RMCOL - 1) +Y 


G LNNO = STARTN(T) + NO 
LMMO STARTM(I) + MO 

See + NO 

ENDM(I) + MO 


LCOL =MAXO STE EEO) 
HCOL =MINO(ENDM(I),HMMO) 
LROW =MAXO STARTN(I), TLNNO) 
HROW =MINO(ENDN(I), HNNO ) 


DOVI33 ROWL=.1 , 3 
RRN = RN + ROW] 
DOtZ5sn, COE = 1 is 
RRM = RM + COL1 
SUM = 0.0 
DO 333 ROW = LROW , HROW 


II=1,K) 


AAN 


ANAAN 


LG 
= 
Z 
=) 
a 
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QAAAaAN 


aAANd 


331 
330 


37 
36 


NNO = ROW - NO 
Ve (NNO ‘LT. LROW) THEN 
RROW = MAXO(ROW,NNO) 
ELSE IF (NNO .GT. HROW) THEN 
RROW = MIN (ROW ,NNO) 
ELSE 


RROW = NNO 
END IF 


DO 433 COL,= LCOL , HGEeL 
MMO) =, COL. = MO 
TE eCMOn it. BC oaye Oe 
CCOL = MAXO (MMO , L) 
ELSE IF (MMOUTGCT. HCOL THEN 
CCOL = MINO (MMO , COL) 


ELSE 
CCOL = MMO 
ENDaGE 
SUM = SUM + FINPUT (ROW, COL,ROW1) * FINPUTSAVEW,CCOL,COL1) 
CONTINUE 
CONTINUE 
R(RRN,RRM) = SUM / ISIZE 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 


THE FOLLOWING FORMAT MUST BE MODIFIED TO USE DIFFERENT FILTER SIZE 
THEN 2 by 2 PIXELS. 


DO’ 3s0 Lb = 71,7 FRG 
WRITE(* , (R(L,J) ,J=1,PKQ) 
WRITE Z 

Soca 12F6.0) 

FORMAT(! ! 

CONTINUE 


RESET THE IDENTITY MATRIX. 
Domeer os = eer 


DO 37 L=1 
cr (jeser a THEN 
MAT The 
IMATRY Gl i= 040 
D IF 
CONTINUE 
CONTINUE 


c 
- RESET SI(J,L) TO HAVE [ I ] IN FIRST PARTITION AND [ O ] IN ALL OTHERS 


AAANAN 


oo 
38 


DO 38 J=1, PK 
DO 39 
IF (J.EQ.L) THEN 
SI(J,L) = 1.0 


SiC, LE) = 700 
IF 


CONTINUE 
CONTINUE 


SOLVE EQUATION [RR] *[B ] 


[ = [ SI - NOTH THAI THE BELOW 
CALL TO IMSL ROUTINE LEQT2F, [ B ] ] 


] 
RETURNS IN [ SI 
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AAAAN OA 


46 
45 


AAAA A 


49 


48 
47 


491 


23 


QAAAA 


Q 


495 


5 


Ui 


IDG = 3 
CALL LEQT2F (R,K,PKQ,PKQ,SI,IDG,WKAREA, IER) 


SOLVE EQUATION [ BOO ] * [ KW ] = [ I ] FOR COVARIANCE 
MATRIX [ KW ] AFTER COLLING IMSL ROUTINE LEQT2F, [ KW ] 
RETURNS IN [ IMATRX ]. 


DO 45 J L K 
DO “46 a 


1 
BOO(J, 35)" = Seles) 
CONTINUE 
CONTINUE 
CALL LEQT2F(B00,K,K,K,IMATRX, IDG,WKAREA, IER) 


SOLVE FOR THE FILTER COEFFICIENTS t ieee) * [ KW), 
WHICH IN THE PROGRAM IS [ A ] = [ SI] * [ IMATRX ] 


DO 47 i eeekO 

Oucces = 1. Kk 
TEMP = 0.0 
pedo Jig =, K 

TEMP = TEMP + SI(J,JJJ) * IMATRK (III, con 

CONTINUE 

A(J,JJ) = TEMP 
CONTINUE 
CONTINUE 


WRITE (*,491) ((TMATRE (J, oye JI=le)), T=) 
FORMAT(' ' , <K>(F6.2,3X)) 
WRITE(*,53) I 
FORMAT (! FILTER COEFFICIENTS :', I2,$) 
WRITE(*,491) ((A(J,3J) , JJ=1,K) , J=1,PKO) 


WRITE OUT MEAN, COVARIANCE MATRIX, AND FILTER COEFFICIENTS TO THE 
USER INPUT FILE. 


OPEN (UNIT=2,FILE= FILTER(T), STATUS='NEW! , CARRIAGECONTROL='LIST', 
FORM='FORMATTED ' ) 


WRITE(2,495) (MEAN (I, J),J=1,K) 

FORMAT (F10. 4 
pre ee 495) a nN JJ=1 ee eee) 
WRITE(2,495 A(J,JJ),JJ=1,K) ,J=1,PKQ) 


CLOSE (UNIT=2) 


CONTINUE 
WRITE (*,55) 

FORMAT(' THE PROGRAM FLTR1 IS OVER',S) 
STOP 
END 


AAANAANANANAAANANANANANANANANAAANAAAANANAAAAAANANAAAA 


On 
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KEKEKKKAKKKARKKKAKKKKAAKAKKRKRRKRKRRRRKRKRRKRKRRKRRKRKRRRRRKRRKRRARRRKRRRRRARKRRRRRKERKAKRARKRRKKI 


PROGRAM SGMT1 


PURPOSE To segment a color image with two textures given 
the image, the image dimensions) thew cer dimen- 
oi 0neaaaa the two sets of filter parameters to 

e used. 


REQUIRED IMSL ROUTINES 
LINV3F, LUDATN, LUELMN 


IMPLEMENTED BY LTJG TIMUR KUPELI Nov 1986 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
KRARKRKKRKKRRKRKRKRKRRKRKRRKRRRRRKRRKRKRRRRRKRRRKARKRRKRKRKRRRARRKRRKRRRKRRKAKRRKRKRRKRERERRERKI 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


RAK VARIABLE DEFINITIONS KKKK 


BINPUT = [{ Input ] Image data in the byte format. | 

ML = { Output The result of ML segmentation in byte formic. 
MAP = { Output _j The result of MAP segmentation in byte format. 
FNAME = [| Input ] Filename of filter parameters set. 

IMAGE = Filename of the image channel. 

Poe 9 = Rows , Columns of filter. : 

N , = Rows , Columns of image. 

K = Number of channel of image. | 

KW = Input |] The covariance matrix... 

MEAN = Input | Mean vector of estimation windows. 

ERROR = Prediction error estimation | 

A = { Input ] Filter coefficients matrix. 

TEXTURE = Zero-mean image data in real*4 format. _ 

Col = Number of removed false points in the first texture. 

Cie = Number of removed false points in the other texture. 

IN , IM = Maximum number of rows and columns in the image. 

Pl , Ql = Maximum number of rows and columns in the filter. 

TMAX = Maximum number of filters for processing. 

IK = Maximum number of channels in image. 


INTEGER IN,N,IM,P1,P,01,0,1TMAX,1T,PO,ROW,COL,1,J,J3J,Jdd , Lee 
LLL,LLLL,COUNT,LI,HI,LJ,HJ,K,PP, OG) TK IL til edly etme 


REAL KW(1:3,1:3,1:2),MEAN(1:3 EE UC einieae ae 


x ERROR(1:128,1:128,1:3,1:2),PML1,PML2, PML LelZ267 Zoo 
* AREA,KS,A(1:3,1:3,1:2,1:2,1:2) ,D1,D2,KW1(1:3,1:3) ,KW2(3,3), 
* EKW1(1,1:3) ,EKW2(1,1:3),PML22,DD2,TEXTUR(1:128,1:128,3,3) 
* WKAREA(6),AA(1:128,1:128) 
CHARACTER*50 IMAGE (1:3),FNAME(1:2) ,MLTEST,MAPTEST 
BYTE BINPUT(1:128,1:128,1:3) ,ML(1:128,1:128) ,MAP(1:128,1:128) 
* MEPCLe1Z5 71 -8Ze) 
IN = 128 
IM = 128 
Pl =4 
1= 4 
MAX = 2 
IK = 3 
GET THE INPUT PARAMETERS OF THE PROGRAM 
WRITE(*,2) IN 
2 FORMAT(' ENTER THE NUMBER OF ROWS IN IMAGE.LIMIT OF',I3,':',$) 
READ(*,3) N 
FORMAT (I3) 


TFC(N.LT.1) sOR. (N2GReIN) Gomees 


A 


52 


ann 


20 


24 
c 


O 


MqQaAMa A A 


Z3 


19 


26 


27 
260 


28 


WRITE(*,5) IM 

FORMAT ( 3) PoeewUnSER OF COLUMNS IN IMAGE.LIMIT OF' ,I3, ,$) 
RE 
PeG@iebt «tj ..OR. (M.GT.IM)) GOTO 4 


WRETE (* 7) 9P1 


FORMAT (! ee THE NUMBER OF ROWS IN FILTER.LIMIT OF',I3, 7S) 

REA 

IF((P.LT.2) .OR. (P.GT.P1)) GOTO 6 

WRITE(*,9) Ql 

+a ENTER THE NUMBER OF COLUMNS IN FILTER.LIMIT OF',I3, Fo) 
R 


3 
MEO, LEC ORE O.GD.Ol)) GOTO 8 


Werte (*~) 1) 7 

FORHAT( ENTER NUMBER Coe texcURboenOePROCESS DIMIT OF" 13,':',$) 
REA 
EOC t.2 )aeOhw (beGT.IMAX)) GO,TO 10 


WRITE(*,13) IK 
FORMAT (' eee Cope UM ee Rh OnmMUMAGESCHANNELS.LIMIT OF',' ',13,$) 
TPR@UK-.LT.1) .OR. (K.GT.IK)) GO TO 12 


GET ALL CHANNELS OF THE IMAGE 


pom f = teak 
WRITE(*,20) I 

FORMAT (' ENTER FILENAME OF THE IMAGE CHANNEL PUHEBR = 12,9) 
READ(*,21) IMAGE(I) 

FORMAT (A50) 


OPEN(UNIT=1, BEE SHAE (1) STATUS='OLD' ,ACCESS='DIRECT' ) 


DO 23 ROW = 
READ (1 ROW) ATEN (RO COLNE) COL = 1 , M*) 
CONTINUE 


CLOSE( UNIT = 1 ) 
CONTINUE 
GET THE FILTER COEFFICIENT,MEAN, AND COVARIANCE MATRICES 


DO 25 I=1,T 
WRITE (*, 26) 

FORMAT('ENTER FILENAME OF FILTER PARAMETERS SET NUMBER',' ',12,$) 
READ(*,21) FNAME(I) 


OPEN (UNIT=2 , FILE=FNAME(1I) ,STATUS='OLD! , FORM='FORMATTED' ) 


DO 260 J = 

READ (2, 27), MEAN( J, i) 
FORMAT(F10.4) 
CONTINUE 


DO 28 J=1,K 
READ(2,27) (KW(J,JJ,1),JJ=1,K) 
CONTINUE 


DO 29 PP=1,P 
DO 30 QQ=1 
DO 31 ROW=1,K 
READ (2, 33) (n (ROW, COL,PP,0Q,I),COL=1,K) 
FORMAT (F10.4) 
CONTINUE 
CONTINUE 
CONTINUE 


5 


amnA AN 


ANAANN 


AAN 


CLOSE (UNIT=2) 


a2 FORMAT (3(F10.4,3X) ) 
Zo CONTINUE 
GET THE OUTPUT FILENAMES OF THE ML AND MAP SEGMENTATION RESULTS. 


aaa ao METEST 
READ(*,220) MAPTEST 
220 FORMAT (A80) 


CONVERT A 2-D BYTE INPUT IMAGE DATA ARRAY IN THE RANGE OF -128 TO 
fae . REPRESENTS APPROPRIATE INTENSITY LEVELS IN THE RANGE OF 
TO : 


N 
DO 37 COL=1,M 
TEMP = BINPUT(ROW,COL,J) 
IF(TEMP.LT.0.0) THEN 
TEMP = TEMP+256 
END IF 


TERT UE Ronee on = TEME = oy 
TEXTUR(ROW,COL,J,2) = TEMP - MEAN(J,2 
37 CONTINUE 
36 CONTINUE 
35 CONTINUE 

CALCULATION OF ERROR ESTIMATE FOR TWO TEATURES 


DO 4091 = tae 
O 41 


or ~LE.0) J=1 
IF(JJJ.LE.0) JJJ=1 


ERROR(L,LL,KK,1)=ERROR(L,LL,KK,1)+A(I1I,JJ,I1I,LLL,1)*TEXTUR(J,JJJ,KK, 1 


46 CONTINUE 
45 CONTINUE 

44 CONTINUE 

43 CONTINUE 

421 CONTINUE 

42 CONTINUE 

41 CONTINUE 

40 CONTINUE 


DO 47 JJ=1,K 
DO 48 LL=1,K 


ae = ee 
KW2(JJ,LL) = KW(JJ,LL,2 
48 CONTINUE 
47 CONTINUE 
Di = 1.0 
CALL _LINV3F (KW1, 6,1,K,K,D1,D2,WKAREA, IER) 
Det Dil 2kAD 2 
LN(1) = avec ene 
Dl = 1.0 


CALL LINV3F(KW2, 6,1,K,K,D1,DD2,WKAREA, IER) 


54 


QAAN 


53 
a2 


54 


iL 


On 


50 


AAA AN 


Det2 = DIle* 24*DDZ 
LN(2) = ALOG(DET2) 


CALCULATION OF MAXIMUM-LIKELIHOOD IMAGE SEGMENTATION 


See N GUNS 3,F ILE=MLTEST ,STATUS='NEW' ,ACCESS='DIRECT', 
ek (IM/4), MAXREC= IN) 


DO 50 ROW =1,N 
BO 5SiecOL = 1,M 
DO S27ie— 1K 


1) =EKW1 (1,1) +ERROR (ROW, Gone AYRES) 
=EKW2(1,1)+ERROR(ROW,COL,J,2)*KW2(J,1 
CONTINUE 

CONTINUE 

PML11 = 0.0 

PML22 = 0.0 

DO 54 L=1,K 

PML11= a ead COLeiy: 2) 


PML22=PML22+EKW2 *ERROR(ROW,COL,L,2 


PMigz2 + ENiCZ 


ML (RO OW, ,COL)=0 
PML(ROW,COL) = PML2 - PML1 
LF(PML1 GT. eee) THEN 
ML(ROW,COL)= - 
DeLr 


MLI(ROW,COL) = ML(ROW,COL) 


= 0:0 
PML1 = PML11 + LN (2) 


CONTINUE 


WRITE(3'ROW) (ML(ROW,COL) , COL=1,M) 
CONTINUE 


CLOSE (UNIT=3) 
MAXIMUM A POSTERIORI IMAGE SEGMENTATION 


OPEN(UNIT= ae eae MAPTES 51a) US= NEW! ,ACCESS="DIRECT', 
RE (IM/4), MAXREC = IN) 


KS = 10.0 

COl = 0 

C10 = 0 

DO 600 II=1, 5 

DO 60 I=1,N 

DOvoL.J = 1, M 

SUM1 = Q.0 
SUMZ = 0.0 
i ees 
pho =i +33 
ie =5 =a 
hoe J 3 


tH 

71 

iy 

Cy 

rm 

O 

Cc 

Cy 
Won ut 
eeteit 


IF(HJ.GT.M) HJ 
PERT Aen nT Fone (ao - Loot 1) 
DO 62 ROW = LI , HI 

DO 63 COL=LJ , Hd 


55 


SUM1 = SUM1 - MLI(ROW,COL) 


63 CONTINUE 
62 CONTINUE 
SUM1 = ae + ME Cie) 
MAP(I,J)= 0 
C 


SUM2 = re 3} - ((KS/AREA)*(2*SUM1-areat1l) ) 
IF(SUM2.LT.0.0) THEN 
MA te =) 
END IF 
MLI(I,J) = MAP(I,J) 
C 
61 CONTINUE 


60 CONTINUE 
600 CONTINUE 


DO 70 IT =1 i N 


? 


DO 71 
1F(ML(T, J) .EQ. 0) T 
IF (MAP (T, rn) -NE. MLCT, J) ) 1GOle= Cole 1 
ee .NE. ML(I,J)) C10=C10+1 
F 
a CONTINUE 
WRITE(4'I) (MAP(I,J),J=1,M) 
70 CONTINUE 


CLOSE (UNIT=4) 
Naeger ,64) cole ,C10 
64 ROBES COl 5x, ClO is) 
O 


END 
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APPENDIX E 


PROGRAMS FOR K-L TRANSFORMATION, AND ONE-CHANNEL SEGMENTATION 


ANAANIAANANANANANANANANNAN 


OQ 


AAAN 


aQ ANN 


100 
101 


102 


(10% 


KREKKKKRKKKRKKRKRKKRKRKRKEKRKRKRRKRKRKRRKRKRKRKRKRKRKRKRKRRKRRRKKRRKRKRKRKRRKRRKRKRRRKRKRKKRKKRRKERKRREKRKEK 


OO OO OE Ob a a OO 


PROGRAM KRLV 


PURPOSE To implement Karhunen-Loeéve transformation from a 2-D 
color image which size is 128 by 128 pixels. 


* 

* 

* 

* 

* 

REQUIRED IMSL ROUTINES 2 
BEGRS;e@ent2s, EHOBKS, EHOUSS, UERTST, USPRKD, UGETIO i 
IMPLEMENTED BY LTJG TIMUR KUPELI Dec 1986 a 
* 


KREKKKRKRKKKEKRKKRKEKRKKRKRKRKRRRKRKRKRKRKRKRRERKRKRKERKRRERKRKRERKRRKRKRKRKRRRKRKRRKEKRKRREKRRERRKRE 


TPEGGR wl Jol kh NL, NZ NN; ROW, COL, ITEMP 0c(128, Zee 
MAXVAL , MINVAL , ‘IDIF, INTVAL, N 


REAL Bids Pleo eh: >, | >3)..F INPUT(1:128,1: eee ese 
Bvs ), WK(3), TEMP ,O(1: :128,1:128,1: pe ‘SLOP 


CHARACTER*50 IMAGE(1:3) , FNAME(1:3) 
BYTE Bleun iel26, 12128) 00(1:128,1:128,1:3) 


TYPE 100 

FORMAT( 'ENTER THE NUMBER OF ROWS, COLUMNS IN THE IMAGE',$) 
READ 101,N 

FORMAT (13) 


Pepe ee? 
FORMAT (" ENTER THE NUMBER OF CHANNELS IN THE IMAGE! AD) 
READ 101,K 


GET THE FILENAMES OF THE RED,GREEN,AND BLUE COMPONENTS 
OF THE COLOR IMAGE. 


NOnl t= lee 
WRITE(*,2) i 
FORMAT(' ENTER THE FILENAME OF THE IMAGE CHANNEL NUMBER ',' ',12,$) 
READ(*,3) IMAGE(I 
WRITE : /3) IMAGE(I) 
WRITE (*.45 
FORMAT (A50) 


CONVERT THE IMAGE BTYE FORMAT TO THE REAL NUMBER FORMAT 
OPEN(UNIT=1, FILE=IMAGE(I), STATUS='OLD!, ACCESS='DIRECT' ) 


DO 5 ROW = 1 
READ (1' ROW) (BINPUT(ROW,COL) , COL=1,N) 
DO 6 COL = 1 
DeMe BINPUT(ROW, gor) 
IF (TEMP. LT. 0.0 an 
TIE 


= TEMP : * 56 
END IF 
FINPUT(ROW,COL,I) = TEMP 
CONTINUE 


CONTINUE 


AAA A NH 


AQAA A AANNA 


23 
22 


21 
20 


59 
30 


Sa 


49 


45 


41 


42 
46 


43 
46 


CLOSE(UNIT = 1) 


CONTINUE 
CALCULATE THE CORRELATION MATRIX 
NN = N * 
DO) 20s eee 
DO: 2s — een 
R(I,J) = 0 
DO 22eNiu = ent 
DOMestice— ee 
R(I,J) = R(I,J) + FINPUT(N1,N2,1I)*FINPUT(N1,N2,J) 
CONTINUE 
CONTINUE 
R(I,J) = R(I,J) / NN 
CONTINUE 
CONTINUE 
WRITE(*,45 
WRITE(*, 30 
WRITE (* , 39 
FORMAT (! ----n-nn nnn nnn nnn nnn nnn nnn n nnn e-e- 
FORMAT sok | THE CORRELATION MATRIX’ ,$) 
WRITES, 31 (RG, I) sae oe) 
FORMAT ( cK F9.2,4X)) 
WRITE (* 


CALCULATE EIGENVALUES AND EIGENVECTORS OF THE 
CORRELATION MATRIX 


JOBN = ll 
CALL EIGRS(R,K,JOBN,D,E,K,WK, IER) 

SORT THE EIGENVALUES IN DECREASING ORDER 
TEMP ey 
aE Des 
Dis TEMG 
DO 49 IT=1, K 

TEMP 


=E 
oe =E 
C3) 2 
CONTINUE 


He 33} 
WRITE(*,45 
FORMAT(' ') 


es. ba 41 
WRITE(*, 39 
FORMAT(' ',5X,!' THE 
DO 46 J=1,K 
WRITE(%, 42) 509) 
FORMAT (5X,F 
CONTINUE 
WRITE(*,39 


EIGENVALUES ',$) 


39) 
FORMAT(! ‘! ' THE EIGENVECTORS ' , $) 
WRITE (* aay CECT. J), J=1,K),I=1,K) 
ei ) 

WRIT 30)" 


USE THE TRANSPOSE OF THE EIGENVECTORS MATRIX TO IMPLEMENT 


38 


QIAN 


QIAN 


63 
62 


KARHUNEN-LOEVE TRANSFORMATION. 


DO 53 J =1 
Q(N1, N2, 1) Q(N1,N2,I) + E(J,1I)*FINPUT(N1,N2,J) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 


D060 I=1 
WRITE(*,61) I 
FORMAT(' ENTER FILENAME OF THE TRANSFORMED IMAGE ',' ',12,8) 
READ(*,3) FNAME(I) 
wae 2) /3) FNAME(I) 
WRITE(*,45) 


CONVERT THE TRANSFORMED IMAGE TO BYTE FORMAT 
DO 62 Nl = 1 N 
DO 63 N2 = 
oc (iy, N2) = JNINT(Q(N1,N2,I)) 
CONTINUE 
CONTINUE 
SCALE THE TRANSFORMED IMAGE TO BE WITHIN DISPLAY RANGE 


ig me 1) THEN 
ek 


MAXVAL = gett: 

MINVAL = 0C(1,1 

DO ROW =1,N 
DO! COL, = 


IF “(Golaon, COL), ST. _UAXVAL) THEN 


Boe le oc abi co COL) ~LT. MINVAL) THEN 
MINVA OC (ROW, COL) 


ENDDO 


INTVAL = MAXVAL - MINVAL 
SLOPE = eae / REAL(INTVAL) 
N 


DO COL 2 =] 
IF (oc (ROW, COL) BQ. "2 MINVAL ) THEN 


OC(R 

Boe IE (Qe(ROW, COL) £9, eS THEN 
OC(ROW, COL) = 

ELSE 


IDIF = QC(ROW,COL) - MINVAL 
QOC(ROW,COL) = INT(SLOPE*IDIF) 


END IF 
ENDDO 
ENDDO 
ELSE 
MINVAL = QC(1,1) 
DO ROW = I h 
DO COL = N 
IF (Q c(RoW, COL) .LT. MINVAL) THEN 
MINVAL '= OC (ROW, COL) 
END IF 
ENDDO 
ENDDO 
DO ROW = N 
DOMcoL y= 1), 


N 
IDIF = QC(ROW,COL) - MINVAL 


a7 


65 
64 


66 
C 


c 
60 


QOC(ROW,COL) = INT(SLOPE*IDIF) 


DO 64 Nl ., N 
DO 65 N2 
ar ENP = QC(N1,N2) 
IF (ITEMP .GT. 134) THEN 
ITEMP = ITEMP - 256 


END IF 
QQ(N1,N2,1) = ITEMP 


"ou 
rap 


CONTINUE 
CONTINUE 


OFEN (UNS ee Bi FNAME (I) ,STATUS='NEW' ,ACCESS='DIRECT', 
(N/4), MAXREC=N ) 


DO 66 Nl =1 

WRITE(2'N1) (Q b (NL, N2,1),N2=1,N) 
CONTINUE 

CLOSE(UNIT= 2) 
CONTINUE 


SLOP 
END 


60 


ANANANNANINAINANNAANAANANAANN 


© 


AQ0AA A AN A 


QAAN 


AAN 


KEKAKKKKKKKKKKKKKRKRKKRKKKKRKKKRKRKRKKRKKRKRKKKKKKKKKKKKRKRKRKKRKRKRKKRKRKRKRRKKKAKKRKKRKRKKKKKKRAKKK 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


10 
la 


2 
13 
14 


15 
16 


PROGRAM FLTR2 


PURPOSE To develop two sets of filter parameters from a 2-D single 
channel image which size is 128 by 128 pixels. These para- 


meters are the mean, the covariance, and the filter 
coefficients. 


REQUIRED IMSL ROUTINES 
LEQT2F, LUDATN, LUELMN, LUREFN 
IMPLEMENTED BY LTJG TIMUR KUPELI Sep 1986 


KRKEKKKKKKKKKKAKKKKKKRKKKKRKRRKRKKKRKAKRKKKKRKRKARRKRRRKAKKKKKKRKKKKKKKKKRKRKKRKKARKAKKKK 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


INTEGER*4 P,Q,ROW,COL,STARTN(2) ,STARTM(2) ,ENDN(2),ENDM(2) , IDGT, 


ROW, RNCOL , RROW, CCOL,X, Ne, 


NO, HMMO , HCOL,HROW, LNNO /LNMO, LCOL,LROW,RN,RM,IER 


REAL*4 SUM, TEMP ,FINPUT(128, 128), MEAN (2) WKAREA(28), 
R(4,4),S 4, Ly, /IMATRX(1,1),B 
KW(1,1),A(4,1),ISIZE 

CHARACTER*50 IMAGE , FILTER(2) 

BYTE BINPUT(128) 

T=2 

GET PROGRAM INPUT PARAMETERS. 

TYPE 

rORHAT(' ENTER NUMBER OF FILTER ROWS FROM 2-4 DESIRED :',$) 
READ 11,P 

FORMAT (I3) 

TYPE 12 

FORMAT(' ENTER NUMBER OF FILTER COLUMNS FROM 2-4 DESIRED :',$) 
READ 11,0 
TYPE 13 

FORMAT(' ENTER NUMBER OF ROWS IN IMAGE ',8) 

READ 11,N 

TYPE 14 

FORMAT(' ENTER NUMBER OF COLUMNS IN IMAGE ',$) 

READ 11,M 

GET FILENAME OF SINGLE CHANNEL IMAGE 

TYPE 15 

FORMAT(' ENTER FILENAME OF IMAGE ',$) 

READ 16,IMAGE 

FORMAT (A50) 

CONVERT THE IMAGE BTYE FORMAT TO REAL NUMBER FORMAT 
OPEN (UNIT=1 1, FILE=IMAGE , STATUS='OLD' ,ACCESS = 'DIRECT') 
DO 1 ROW =1, N 

READ | (1' ROW) _(BINPUT (COL), ,COL=1,M) 
TEMP = ‘BINPUT (COL) 
IF (TEMP.LT.0.0) THEN 
TEMP = TEMP + 256 


END IF 
FINPUT(ROW,COL) = TEMP 


61 


AAAN 


aan 


aAaAN 


AAA 


18 
17 


20 


ZL 


22 


23 


24 


26 
25 


Zi 


20 
28 


291 


CONTINUE 


CONTINUE 
CLOSE (UNIT=1) 


GET 


'T' AREAS FOR WHICH FILTERS ARE DESIRED AND OUTPUT FILENAME FOR 


EACH AREA'S FILTER COEFFICIENTS AND COVARIANCE MATRIX. 


DO 


1S = 155 
WRITE(*,20) I 
FORMAT(' ENTER UPPER-LEFT ROW FOR AREA ',12,':',$) 


READ(*,11) STARTN(I) 
WRITE(*,21) I 


FORMAT (' ENTER UPPER-LEFT COLUMN FOR AREA! , 12, ' 2555) 


READ(*,11) STARTM(I) 
WRITE (4322) ee 


FORMAT (' ENTER NCI) RIGHT ROW FOR AREA' ,I2,':',$) 


READ(*,11) ENDN(I 
WRITE(*,23) I 


FORMAT(' ENTER LOWER-RIGTH COLUMN FOR AREA! ,I2,':',$) 


READ(*,11) ENDH(T) 
WRITE(*,24) I 


FORMAT (' ENTER OUTPUT FILE-SPEC FOR FILTER! p25 


READ(*,16) FILTER(I) 


FIND THE MEAN VECTOR -OF ESTIMATION WINDOW AREA OF IMAGE 
1S1ZE : (ENDN (I) - STARTN(I) + 1)*(ENDM(I) - STARTM(I) + 1) 
DO 25. L= STARIN(T) , ENDN(I | 
DO 26 


) 
J = STARTM(I), ENDU(T) 
SUM = SUM + FINPUT(L,J) 
CONTINUE 


CONTINUE 

MEAN (I) = SUM / ISIZE 

WRITE eS poe aS 

FORMAT MEAN! , (122) eno coe 


CORRECT THE IMAGE TO BE ZERO MEAN 


DOVZ6ee Ge Se ENDN (I 
DO 29 


) 
J = TARTM (I) , ENDM(TI) 
FINDUT(L, J) = FINPUT(L, J) - MEAN(T) 
CONTINUE 


CONTINUE 
DETERMINE CORRELATION MATRIX 


WRITE (* ,291) I 
FORMAT 


=o 


PQ Q 
DO 30 RNROW = 


CORRELATION.MATRIX'! ,I2,$) 


P 
X = Q * (RNROW-1) 
DO 31 RNCOL = 1,P 
NO_= RNROW - RNCOL 
a ous (RNCOL-1 ) 
DO 32 RMROW = 1,0 
RN = X + RMROW 
DO 33 RMCOL = 1,9 
MO = RMROW - RMCOL 


RM = Y + RMCOL 

LNNO = SR + NQ 
LMMO = STARTM(I) + MO 
HNNO = Se + NOQ 
HMMO = ENDM(1I) + MO 


ON 
i] 


C 
55 
34 
Cc 
33 
32 
a 
30 
c 
C 
C 
41 
c 
e 
c 
c 
37 
25 
36 
c 
C 
c 
C 
360 
Ec 
S 
c 
€ 
e 
Cc 
C 
c 
43 
C 
C 
C 
e 


LCOL =MAXO(STARTM(I) ,LMMO) 
HCOL =MINO(ENDM(1) ,HMMO) 
LROW =MAXO(STARTN(I),LNNO) 
HROW =MINO(ENDN(1I),HNNO) 


SUM = 0.0 
DO 34 ROW = ee HROW 


NN ROW - NO 
DOrss COL = LCOL, HCOL 


MMO = COL - MO 
SUM = SUM + FINPUT(ROW,COL)*FINPUT(NNO,MMO) 
CONTINUE 
CONTINUE 
R(RN,RM) = SUM / ISIZE 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 
RESET SI(J,1) TO HAVE 1 IN THE FIRST ROW AND O IN ALL OTHERS 
DO 41 


a= eo 
IF (J.EQ.1) THEN 
eno 
IMATRX(J,1) 
Gu(ael = 0.0 


lo 


eo 


CONTINUE 


THE FORMATS BELOW MUST BE MODIFIED TO USE DIFFERENT FILTER 
SIZE THAN 2 by 2 PIXELS. 


IDGT = 1 
DO 36 i 4h0 
WRITE a (R(K,J),J=1,PQ) 
WRITE 38 
FORMAT ,4F9.2, $) 
FORMAT (' ES) 
CONTINUE 
SOLVE EQUATION [RJ] *[{BJ=[ SI ] 
. CALL LEQT2F (R,IDGT,PQ,PQ,SI,3,WKAREA, IER) 


WRITE(*, 360) (SI(J,1),J=1,PQ) 
FORMAT (F10.2) 


SOLVE EQUATION B00 * KW =I 
BOO = SI(1,1) 
GET THE FILTER COEFFICIENTS USING THE EQUATION 
[A ]=[B] * kW 


KW(1, 1) = IMATRX(1,1) / BOO 
ee ,*) ‘COVARIANCE! 
WRITE(*,360) KW(1,1) 
DO 43 ROW = 1,P9 

A(ROW,1) = £1 (ROW, 1) * KW(1,1) 
CONTINUE 


WRITE OUT MEAN, COVARIANCE, AND FILTER COEFFICIENTS TO THE 
Use oN UT FILE . 


OPEN (UNIT=2,FILE=FILTER(1I) ,STATUS='NEW' , CARRIAGECONTROL='LIST'! 
FORM=!FORMATTED' ) 


50 
500 

aL 

52 


53 
54 





WRITE(2,50) KW(1,1) 
FORMAT( F10.5 
WRITE 2,500) MEAN (I) 
FORMAT( F10.5) 
DO 51 J = 1,PQ 
WRITE(2,52) A(J,1) 
CONTINUE 


FORMAT (F10.5) 
CLOSE (UNIT=2) ° 


WRITE(*,53) I 

FORMAT(!' FILTER COEFFICIENTS :',12,$) 
WRITE(*,54) (A(J,1),J=1,P0Q) 
FORMAT(F10.5) 

CONTINUE 

STOP 

END 
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ANAANANANAANNANANAAARAN 


OQ 


AAN 


20 
2a 


KKAKRKKRKKRKRRKRKRKKKRRKARKRKRRKKRKRRRRKRKRRRRKRRRRKRRRRKRRKAKRKRRRKRKRRRRKRRKRKRKRRRKRKRRRKKKKRKK 


PS i a 


PROGRAM SGMT2 

PURPOSE To segment a single image with two textures given the 
image, the image dimensions, the filter dimensions, 
and two sets of filter parameters to be used. 

REQUIRED IMSL ROUTINES 
NONE 


PIP EEMENTED BY LIJG TIMUR KUPELI sep 1986 


REKKKKKKKRKKRKKKRKKRRKKRKAKRKRKRKRRRRRKRERKRKRRKKRRKRRKRKRRRRRRERRRRRRKRRRKRRKRKRRKKR 


INTEGER IN,N,IM,P1,P,01,0,TMAX,T,PQ,ROW,COL,I,J,JJ,JJJ,L,LL, 
UEUSEDLE COUNT, LL,HI,LJ,HJ 


REAL KW(2),MEAN(2),AA(4,1,2),TEMP,SUM1 ,SUM2,TEXTUR(128,128,2), 
LN(2),ERROR(128,128,2),PML1,PML2,PML(128,128), 
AREA,KS ,A(2,2) 

CHARACTER*50 IMAGE , FNAME ,MLTEST,MAPTEST 


BYTE PINEUi2ce2s ME 128-128) MAP i128, 128) ,MLI(128,128) 


IN = 128 
IM = 128 
Pl = 4 
1= 4 
MAX = 2 


GET THE INPUT PARAMETERS OF THE PROGRAM 


WRITE(*,2) IN 

FORMAT(' ENTER THE NUMBER OF ROWS IN IMAGE.LIMIT OF',I3,':',$) 
READ(*,3) N 

FORMAT (I3) 
IF((N.LT.1) .OR. (N.GT.IN)) GOTO 1 


WRITE(*,5) IM 
FORMAT ( Bate THE NUMBER OF COLUMNS IN IMAGE.LIMIT OF',I3,':',$) 
[PCG LE1) OR. (M-GT.IM)) GoTo 4 


WRITE(*,7) Pl 


et ENTER THE NUMBER OF ROWS IN FILTER.LIMIT OF',I3,':',$) 
Ea aid ORR AGI PL)) GOTO 6 

WRITE(*,9) Q1 

le 1 EN ER THE NUMBER OF COLUMNS IN FILTER.LIMIT OF',I3,':!',$) 


IF((Q.LT.2) .OR. (Q.GT.Q1)) GOTO 8 


WRITE(*,11) TMAX 
FORMAT, ENTER NUMBER OF TEXTURES TO PROCESS.LIMIT OF',1I3,':',$) 
Pecirie2) Oen Cl.GT.TMnx)) GO rOnLe 


GET THE SINGLE-CHANNEL IMAGE 
PO = P * 
WRITE (*, 20) 
FORMAT(' ENTER FILENAME OF IMAGE ',$) 


READ(*,21) IMAGE 
FORMAT (A50) 
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* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


QO 


AAA A 


Ze 


24 


240 


G 
242 
c 


AAAA 


241 
ZS 


220 


42 
41 


46 


OPEN(UNIT=1, FILE=IMAGE , STATUS='OLD' , ACCESS='DIRECT' ) 


DO 22 ROW = 
READ (1' ROW) (BINPUT(ROW,COL),COL = 1 , M ) 
CONTINUE 


CLOSE( UNIT = 1 ) 
GET THE FILTER COEFFICIENTS,MEANS, AND COVARIANCES 


pO 23 I=1,T 

WRITE(*,24) I 

FORMAT(' ENTER THE FILENAME OF FILTER PARAMETERS SET NUMBER! ,12,':',$) 
READ(*,21) FNAME 


OPEN (UNIT=2, FILE=FNAME , STATUS='OLD' ) 


READ(2, ne KW(TI) 
FORMAT (F1 
READ(2, 240) MEAN(TI) 
DO 242 K lee ° 
READ (2, 340) AA KK,1,1) 


CONTINUE 
CLOSE(UNIT = 2) 


WRITE(%, cae KW(I),MEAN(I),AA(1,1 it), AA(2,1,1),AA(3,1,1I),AA(4,1,1) 
FORMAT (6F10 
CONTINUE 
eee! ae HEPES. 
READ(*,220) MAPTEST 
FORMAT (A80) 


CONVERT A 2-D BYTE INPUT IMAGE DATA IN THE RANGE OF -128 TO lZ7enaee 
REPRESENTS APPROPRIATE INTENSITY LEVELS IN THE RANGE OF O TO 255 


DO 30 ee 1 N 

© ei ecommam! M 

TEMP = BINPUT(ROW, COL) 

IF (TEMP. Li seo TEMP = TEMP+256 

TEXTUR (ROW, ,COL, z = TEMP = eats 

TEXTUR(ROW,COL,2) = TEMP - MEAN(Z 
CONTINUE 

CONTINUE 


CALCULATOIN OF ERROR ESTIMATE FOR TWO TEXTURES 
DO 40 I=1,T 
J = 


Dosa) onan 
DO 42 JJJ = fo 
ACIS, 333) = AMG oe lege) 
je] 


CONTINUE 
CONTINUE 


DO 43 L=l N 
DO 44 LL = 1 M 
ERROR(L, eee I) = 0.0 
DO 45 LLL =1 , oe 
J = LeeLee 
DO 46 LLLL = 1 , Q 
JJ = LL - LLLL 
IF((J.GT.0) .AND. (JJ.GI20))) 0sEh 
ERROR(L,LL,I) “ERE LE L)+ A(LLL, LLLL)*TEXTUR(J, Jd 7k) 


CONTINUE 


66 


45 CONTINUE 
44 CONTINUE 
43 CONTINUE 


LN(I) = ALOG(KW(I)) 
40 CONTINUE 


c 
C CALCULATION OF MAXIMUM LIKELIHOOD IMAGE SEGMENTATION 
c 


OPEN(UNIT=3 , FILE=MLTEST , STATUS='NEW' , ACCESS='DIRECT', 
RECL=(IM/4) ,MAXREC=IN) 


ad 
nee 
Dersone — 1 
Horst 3 = 1. M 
PML1 = 0.0 
PML2 = 0.0 
PML1 epee a hyena 
PML2 ERROR(1I,J,LL)**2)/KW(L 
PML(I,J) = PML2 - PML1 
IF(PML1 .GT. PML2) THEN 
ML(I,J)= -1 
END IF 
Mend, J) = NL(I, 3) 
51 CONTINUE 


WRITE(3'I) (ML(I,J) ,J=1,M) 
50 CONTINUE 


= 


L) 


On 


CLOSE (UNIT=3) 
MAXIMUM A POSTERIORI IMAGE SEGMENTATION 


OPEN (UNIT=4 , FILE=MAPTEST , STATUS='NEW! , ACCESS='DIRECT', 
RECL= (IM/4), MAXREC = IN) 


KS = 100 

DO 600 II =1, 5 
DO 60 I =1 ,N 

DO 61 J = eel 

S UML 

SUM2 


OQ g0NA aN 


tr 
Cc, 
tt it tt tl 
CCH Il Il 
t 
+1 +100 
WWWWOO 


PE GIVEO.O) Lil 
Peal Gl .N) HE 
IF(LJ.EO.0) LJ 
IF(HJ.GI.M) HJ 


AREA = (HI - LI +1) * (HJ - LJ +1) 


DO 62 ROW = LI , HI 
DO 63 COL=LJ , HJ 
IF((ROW.NE.I) .OR. (COL.NE.J)) THEN 
SUM1 = SUM1 - MLI(ROW,COL) 
END IF 
63 CONTINUE 
62 CONTINUE 
MAP(I,J)=0 


SUM2 = PML(I,J) - 
IF(SUM2.LT.0) THEN 
MAP(I,J) = -1 

ND IF 


Te | | | | 
<r ae 


((KS/AREA) * (2*SUM1L-AREA+1 ) ) 


67 


61 
600 
Cc 


70 


MLI(I,J) = MAP(I,J) 
CONTINUE 
CONTINUE 
CONTINUE 


BO’ 70 I = 1, 7 N 
WRITE(4'I) (MAP(I,J),J=1,M) 
CONTINUE 
STOP 


END 
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